################################################################################
#Read Data																						     
################################################################################
rslt.file = paste(c(rslt.path,"r-mcb-estimated.txt"),collapse = "")

mcb.estimated = array(dim = c(4, 99,3))
for(e in 1:4){
	lower = 99*(e - 1) + 1; upper = 99*e
   mcb.estimated[e,,] = as.matrix(read.table(rslt.file, na.strings = ".")[lower:upper,])
}
mcb.estimated = mcb.estimated/4
################################################################################
#Step 1) Export Results   		                                                  
################################################################################
file.list     = c("mbt","mct","mst")
estimand.list = c(expression(paste(B[I]-MTE,"(",x[I],",",u[S],")")),
                  expression(paste(C[I]-MTE,"(",z,",",u[S],")")),
                  expression(paste(S[I]-MTE,"(",x[I],",",z,",",u[S],")")))
legend.pos  = c("topright","topleft","topright")
par.pch 		= c(15, 25, 22, 10)
xlim.par 	= c( 0.0,1.0)

x.idx_single = seq(1,99,by = 2); u.seq_single = seq(0.01, 0.99, by = 0.01)
x.idx_double = seq(1,99,by = 2); u.seq_double = seq(0.01, 0.99, by = 0.02)
################################################################################
#Step 1b): Export All Results Single Graphs                                    
################################################################################
ylim.par = c(floor(min(mcb.estimated)), ceiling(max(mcb.estimated)))
pdf(file= paste(c(out.path,"g-all-result.pdf"), collapse =""), width = 7,  height = 5)
#-------------------------------------------------------------------------------
#parameter for the graph
#-------------------------------------------------------------------------------
par(mar = c(5,4,4,0), mgp = c(3,1,0))

for(i in 1:3){
   if(i == 1){
      plot(u.seq_double, mcb.estimated[(i + 1),x.idx_double,1], main = "", 
        pch = par.pch[1],xlab = "", ylab= "", col = "black", lwd = 1,
        yaxs = "i", ylim = ylim.par, xlim = xlim.par, xaxs = "i",
        bty = "l");
      lines(u.seq_single,  mcb.estimated[(i + 1),,1])
  
      title(xlab=expression(u[S]), ylab = "Marginal Effects of Treatment")
      legend("bottom", legend = estimand.list[1:3], bty = "n", y.intersp = 1.5,
        lty = c(1,1,1)  , pch = par.pch[1:3], col = "black", cex = 1.0)
   }  else {
      par(new = TRUE)
      plot(u.seq_double,  mcb.estimated[(i + 1),x.idx_double,1], main = "",
        pch = par.pch[i],xlab = "", ylab= "", col = "black", lwd = 1,
        yaxs = "i", ylim =  ylim.par, xlim = xlim.par, xaxs = "i",
        axes = FALSE, bty = "l");
      lines(u.seq_single,  mcb.estimated[(i + 1),,1])
   } 
}
dev.off()
################################################################################
#Step 1c): Export Results Single Graphs                                       
################################################################################
for(e in 1:3){
	ylim.par = c(floor(1.05*min(mcb.estimated[(e + 1),,2])), ceiling(1.05*max(mcb.estimated[(e + 1),,3])))
	pdf(file=paste(c(out.path,"g-",file.list[e],"-result.pdf"),collapse=""),
		 width = 7, height = 5)
   plot(u.seq_double, mcb.estimated[(e + 1),x.idx_double,1], main = "",
        pch = par.pch[1],xlab = "", ylab= "", col = "black", lwd = 1, yaxs = "i",
        ylim = ylim.par, xlim = xlim.par, xaxs = "i", bty = "l")
   lines(u.seq_single, mcb.estimated[(e + 1),,1])
   lines(u.seq_single, mcb.estimated[(e + 1),,2], lty = 2, lwd = 1)
   lines(u.seq_single, mcb.estimated[(e + 1),,3], lty = 2, lwd = 1)
  
   title(xlab=expression(u[S]), ylab = "Marginal Effects of Treatment")
   legend(legend.pos[e], legend = c(estimand.list[e],"90% Confidence Interval"),
          bty = "n", y.intersp = 1.5,  lty = c(1,2)  , pch = c(15, NA),
          col = "black", cex = 1.0)
   dev.off()
}
################################################################################
#Step 1c): marginal support                                      
################################################################################
max.count = 60
rslt.file = paste(c(rslt.path,"r-propensity-scores.txt"),collapse = "")

P = as.matrix(read.table(rslt.file, na.strings = "."))[,1]
D = as.matrix(read.table(rslt.file, na.strings = "."))[,2]

pdf(file = paste(c(out.path,"g-marginal-support.pdf"), collapse =""), width = 7,height = 5)
#-------------------------------------------------------------------------------
#parameter for the graph
#-------------------------------------------------------------------------------
par(mar = c(5,4,4,0), mgp = c(3,1,0))
   marginal.support = list(P[(D == 1)] ,P[(D == 0)])
   multhist(marginal.support,breaks = seq(0.00,1.00, by = 0.02), ylab = "Counts",
            xlab = "Propensity Score",ylim = c(0,max.count)) 
   legend("top", legend = c("Treated","Untreated")
          , bty = "n", y.intersp = 1.5, pch = c(15,15)
          , col = c("black", "grey"), cex = 1.3)
dev.off()
################################################################################
################################################################################






















